######################################################################
######################################################################
# Replication Code for: 
# Woldense and Kroeger. Elite Change without Regime Change:
# Authoritarian Persistence in Africa and the End of the 
# Cold War. American Political Science Review. 
######################################################################
######################################################################

#The code below replicates the analyses, figures, and tables in the 
#main text and supplemental materials. 

#Windows users must first install Rtools
#See: https://cran.r-project.org/bin/windows/Rtools/rtools42/rtools.html

#Install devtools package to install previous versions of Zelig and stargazer
install.packages('devtools')
library(devtools)


#install last version of Zelig (before it was
#removed from CRAN) and its dependencies
install_version('Zelig', version = '5.1.7',
                repos = "http://cran.us.r-project.org",
                dependencies = TRUE)

#install pacman to check if other needed packages
#are installed, install any missing packages, 
#and load all packages
install.packages('pacman')

pacman::p_load(multiwayvcov, lmtest, xtable, 
               margins, survival, car, ggpubr, 
               zoo, tidyverse)

#install previous version of stargazer 
#current version is not compatible with the code
#below
install_version('stargazer', version = '5.2.2', 
                repos = "http://cran.us.r-project.org")

#load stargazer and Zelig
library(stargazer)
library(Zelig)


#load full dataset
df <- read.csv("Woldense_Kroeger_APSR.csv")


#Main GWF sub-Saharan Data Sample ------------------------------------------------------------

#list of countries to omit for main sub-Saharan analyses, includes 
#south africa under apartheid 
nafrica_sa <- c('Morocco','Algeria','Tunisia','Libya','Egypt','South Africa')

gwf_ssa <- df %>%
  filter(!is.na(gwf_regimetype)) %>%
  filter(!(country_name %in% nafrica_sa))


# Figure 1 ----------------------------------------------------------------
#Cabinet Change in Gabon, 1968-2008

gwf_ssa %>%
  filter(cowcode==481, 
         #select years where Bongo led the cabinet observed by WhoGov
         year > 1967 & year < 2009) %>%
  arrange(year) %>%
  ggplot(., aes(x = year, y = minister_tenure, color = as.factor(spell_end))) +
  geom_point(position = position_jitter(.1), alpha = .4) +
  scale_color_manual(values = c('gray60', 'red'),
                     name = 'Minister Status:',
                     labels = c('Survived',
                                'Dismissed')) +
  labs(x = 'Year', y = 'Minister Tenure') +
  scale_x_continuous(breaks = c(1970, 1980, 1990, 2000, 2008)) +
  theme_pubclean() +
  theme(legend.key = element_rect(fill = NA))



# Figure 2 ----------------------------------------------------------------
#Minister exits across sub-Saharan autocracies, 1966-2010

#Yearly percent of ministers exiting across all countries in gwf_ssa sample
minister_turnover <- gwf_ssa %>%
  filter(!is.na(gwf_regimetype)) %>%
  filter(!(country_name %in% nafrica_sa)) %>%
  arrange(cowcode, year) %>%
  group_by(year) %>%
  summarise(minister_exits = sum(spell_end), 
            total_ministers = n(), 
            exit_percent = 100*(minister_exits/total_ministers)) %>%
  mutate(type = 'Minister') %>%
  dplyr::select(year, exit_percent, type)


#Minister turnover across SSA, yearly mean and 5 year moving average

ma5 <- data.frame(exit_ma = rollmean(minister_turnover$exit_percent, k = 5),
                  year = 1968:2008)

minister_turnover_ma <- left_join(minister_turnover, ma5) %>%
  pivot_longer(cols = contains('exit')) %>%
  mutate(name = factor(name, levels = c('exit_percent', 'exit_ma')))

ggplot(minister_turnover_ma, aes(x = year, y = value, 
                                 color = as.factor(name), 
                                 linetype = as.factor(name))) +
  geom_line() +
  labs(x = 'Year', y = '% Exiting Cabinet') +
  scale_x_continuous(breaks = c(1966, seq(1970, 2010, 5))) +
  scale_color_manual(values = c('gray75', 'black'),
                     labels = c('Yearly data', '5-year moving average')) +
  scale_linetype_manual(values = c(2, 1),
                        labels = c('Yearly data', '5-year moving average')) +
  theme_pubclean() +
  theme(panel.background = element_rect(fill = 'white'),
        legend.key = element_rect(fill = 'white'),
        legend.title = element_blank(),
        axis.text.x = element_text(angle = 0),
        panel.grid.major.y = element_blank()) 


# Figure 3 -----------------------------------------------------------
#Average marginal effects of Cold War end window variables


#Figure 3 models and tables of results

#list of cwe variables (5-years - 1-year)

cwe_options <- list('CWE_5', 'CWE_4', 'CWE_3', 'CWE_2',
                    'CWE_1')


#control options

none <- ''


baseline <- '+ ln_oilgas_rents + 
                    pers_2pl +
                    lcgdppc + 
                    cabinet_size + 
                    left_truncated'

tenure <- '+ minister_tenure +
                    I(minister_tenure^2) +
                    I(minister_tenure^3) +
                    ln_oilgas_rents + 
                    pers_2pl +
                    leader_tenure +
                    I(leader_tenure^2) +
                    I(leader_tenure^3) +
                    lcgdppc + 
                    cabinet_size + 
                    left_truncated'


events <- '+ minister_tenure +
                  I(minister_tenure^2) +
                  I(minister_tenure^3) +
                  ln_oilgas_rents + 
                  pers_2pl +
                  leader_exit + 
                  leader_tenure +
                  I(leader_tenure^2) +
                  I(leader_tenure^3) +
                  election_corrected + 
                  coup_corrected + 
                  lcgdppc + 
                  cabinet_size + 
                  left_truncated'

events_no_oil <- '+ minister_tenure +
                  I(minister_tenure^2) +
                  I(minister_tenure^3) +
                  pers_2pl +
                  leader_exit + 
                  leader_tenure +
                  I(leader_tenure^2) +
                  I(leader_tenure^3) +
                  election_corrected + 
                  coup_corrected + 
                  lcgdppc + 
                  cabinet_size + 
                  left_truncated'

aid_dem <- '+ minister_tenure +
                  I(minister_tenure^2) +
                  I(minister_tenure^3) +
                  ln_oilgas_rents +
                  ln_total_aid_AD +
                  v2x_polyarchy +
                  pers_2pl +
                  leader_exit + 
                  leader_tenure +
                  I(leader_tenure^2) +
                  I(leader_tenure^3) +
                  election_corrected + 
                  coup_corrected + 
                  lcgdppc + 
                  cabinet_size + 
                  left_truncated'


model_specs <- list(none, baseline, tenure, events)

model_specs_all <- list(none, baseline, tenure, events, aid_dem)

#Pooled logit models used to calculate AMEs in Figure 3
pooled_logit <- function(cwe_window){
  
  lapply(model_specs, function(z){
    
    model <- glm(as.formula(paste('spell_end ~', 
                                  cwe_window, z)), 
                 data = gwf_ssa, 
                 family = binomial)
    
    
    SEs <- sqrt(diag(cluster.vcov(model, gwf_ssa$leaderID)))
    
    return(list(model, SEs))
    
  })}


#5-year window models 

cwe5 <- pooled_logit('CWE_5')

#Table of results for cwe 5-year models above
#stargazer(cwe5[[1]][[1]],
#          cwe5[[2]][[1]],
#          cwe5[[3]][[1]],
#          cwe5[[4]][[1]], 
#          se = list(cwe5[[1]][[2]],
#                    cwe5[[2]][[2]],
#                    cwe5[[3]][[2]],
#                    cwe5[[4]][[2]]),
#          table.placement = '!hbt',
#          dep.var.labels = 'Minister Exit',
#          covariate.labels = c('CWE 5-year', 'Minister Tenure',
#                               'Minister Tenure$^{2}$', 
#                               'Minister Tenure$^{3}$',
#                               'Ln(Oil/Gas Rents)',
#                               'Ln(Total Aid)', 
#                               'VDem Polyarchy',
#                               'Personalism', 'Leader Exit',
#                               'Leader Tenure',
#                               'Leader Tenure$^{2}$', 
#                               'Leader Tenure$^{3}$',
#                               'Election',
#                               'Failed Coup',
#                               'Ln(GDP per Capita)', 
#                               'Cabinet Size', 'Left Truncated'),
#          omit.stat = c('ll', 'aic'),
#          digits = 2,
#          no.space = T,
#          title = 'Pooled Logit Models of Minister Survival with 5-year CWE Window', 
#          notes = 'Standard errors clustered by leader.', 
#          star.cutoffs = c(0.05, 0.01, 0.001))

#4-year window models

cwe4 <- pooled_logit('CWE_4')

#Table of results for cwe 4-year models above
#stargazer(cwe4[[1]][[1]],
#          cwe4[[2]][[1]],
#          cwe4[[3]][[1]],
#          cwe4[[4]][[1]], 
#          se = list(cwe4[[1]][[2]],
#                    cwe4[[2]][[2]],
#                    cwe4[[3]][[2]],
#                    cwe4[[4]][[2]]),
#          table.placement = '!hbt',
#          dep.var.labels = 'Minister Exit',
#          covariate.labels = c('CWE 4-year', 'Minister Tenure',
#                               'Minister Tenure$^{2}$', 
#                               'Minister Tenure$^{3}$',
#                               'Ln(Oil/Gas Rents)',
#                               'Ln(Total Aid)', 
#                               'VDem Polyarchy',
#                               'Personalism', 'Leader Exit',
#                               'Leader Tenure',
#                               'Leader Tenure$^{2}$', 
#                               'Leader Tenure$^{3}$',
#                               'Election',
#                               'Failed Coup',
#                               'Ln(GDP per Capita)', 
#                               'Cabinet Size', 'Left Truncated'),
#          omit.stat = c('ll', 'aic'),
#          digits = 2,
#          no.space = T,
#          title = 'Pooled Logit Models of Minister Survival with 4-year CWE Window', 
#          notes = 'Standard errors clustered by leader.', 
#          star.cutoffs = c(0.05, 0.01, 0.001))



#3-year window models and table

cwe3 <- pooled_logit('CWE_3')

#Table of results for cwe 3-year models above
#stargazer(cwe3[[1]][[1]],
#          cwe3[[2]][[1]],
#          cwe3[[3]][[1]],
#          cwe3[[4]][[1]], 
#          se = list(cwe3[[1]][[2]],
#                    cwe3[[2]][[2]],
#                    cwe3[[3]][[2]],
#                    cwe3[[4]][[2]]),
#          table.placement = '!hbt',
#          dep.var.labels = 'Minister Exit',
#          covariate.labels = c('CWE 3-year', 'Minister Tenure',
#                               'Minister Tenure$^{2}$', 
#                               'Minister Tenure$^{3}$',
#                               'Ln(Oil/Gas Rents)',
#                               'Ln(Total Aid)', 
#                               'VDem Polyarchy',
#                               'Personalism', 'Leader Exit',
#                               'Leader Tenure',
#                               'Leader Tenure$^{2}$', 
#                               'Leader Tenure$^{3}$',
#                               'Election',
#                               'Failed Coup',
#                               'Ln(GDP per Capita)', 
#                               'Cabinet Size', 'Left Truncated'),
#          omit.stat = c('ll', 'aic'),
#          digits = 2,
#          no.space = T,
#          title = 'Pooled Logit Models of Minister Survival with 3-year CWE Window', 
#          notes = 'Standard errors clustered by leader.', 
#          star.cutoffs = c(0.05, 0.01, 0.001))


#2-year window models and table

cwe2 <- pooled_logit('CWE_2')

#Table of results for cwe 2-year models above
#stargazer(cwe2[[1]][[1]],
#          cwe2[[2]][[1]],
#          cwe2[[3]][[1]],
#          cwe2[[4]][[1]], 
#          se = list(cwe2[[1]][[2]],
#                    cwe2[[2]][[2]],
#                    cwe2[[3]][[2]],
#                    cwe2[[4]][[2]]),
#          table.placement = '!hbt',
#          dep.var.labels = 'Minister Exit',
#          covariate.labels = c('CWE 2-year', 'Minister Tenure',
#                               'Minister Tenure$^{2}$', 
#                               'Minister Tenure$^{3}$',
#                               'Ln(Oil/Gas Rents)',
#                               'Ln(Total Aid)', 
#                               'VDem Polyarchy',
#                               'Personalism', 'Leader Exit',
#                               'Leader Tenure',
#                               'Leader Tenure$^{2}$', 
#                               'Leader Tenure$^{3}$',
#                               'Election',
#                               'Failed Coup',
#                               'Ln(GDP per Capita)', 
#                               'Cabinet Size', 'Left Truncated'),
#          omit.stat = c('ll', 'aic'),
#          digits = 2,
#          no.space = T,
#          title = 'Pooled Logit Models of Minister Survival with 2-year CWE Window', 
#          notes = 'Standard errors clustered by leader.', 
#          star.cutoffs = c(0.05, 0.01, 0.001))


#1-year window models and table

cwe1 <- pooled_logit('CWE_1')

#Table of results for cwe 1-year models above
#stargazer(cwe1[[1]][[1]],
#          cwe1[[2]][[1]],
#          cwe1[[3]][[1]],
#          cwe1[[4]][[1]], 
#          se = list(cwe1[[1]][[2]],
#                    cwe1[[2]][[2]],
#                    cwe1[[3]][[2]],
#                    cwe1[[4]][[2]]),
#          table.placement = '!hbt',
#          dep.var.labels = 'Minister Exit',
#          covariate.labels = c('CWE 2-year', 'Minister Tenure',
#                               'Minister Tenure$^{2}$', 
#                               'Minister Tenure$^{3}$',
#                               'Ln(Oil/Gas Rents)',
#                               'Ln(Total Aid)', 
#                               'VDem Polyarchy',
#                               'Personalism', 'Leader Exit',
#                               'Leader Tenure',
#                               'Leader Tenure$^{2}$', 
#                               'Leader Tenure$^{3}$',
#                               'Election',
#                               'Failed Coup',
#                               'Ln(GDP per Capita)', 
#                               'Cabinet Size', 'Left Truncated'),
#          omit.stat = c('ll', 'aic'),
#          digits = 2,
#          no.space = T,
#          title = 'Pooled Logit Models of Minister Survival with 1-year CWE Window', 
#          notes = 'Standard errors clustered by leader.', 
#          star.cutoffs = c(0.05, 0.01, 0.001))



#function to calculate average marginal effects
#from pooled logit models with
#specifications varying the CWE window (5yr-1yr) 
#and the control variable sets

fig3 <- lapply(cwe_options, function(x){
  
  lapply(model_specs, function(y){
    model <-  glm(as.formula(paste('spell_end ~', 
                                   x, y)), 
                  data = gwf_ssa, 
                  family = binomial)
    
    clustered_vcov <- cluster.vcov(model, gwf_ssa$leaderID)
    
    #marginal effects 
    data.frame(summary(margins(model, vcov = clustered_vcov)))
  })
})


#Figure 3 plot 
bind_rows(fig3) %>%
  filter(grepl("CWE", factor)==T) %>%
  mutate(specification = rep(c('No Controls', '+ Baseline', 
                               '+ Tenure', '+ Events'), 5), 
         specification = factor(specification, 
                                levels = c('No Controls',
                                           '+ Baseline',
                                           '+ Tenure',
                                           '+ Events')),
         window_size = case_when(factor=='CWE_5' ~ 1, 
                                 factor=='CWE_4' ~ 2, 
                                 factor=='CWE_3' ~ 3, 
                                 factor=='CWE_2' ~ 4, 
                                 factor=='CWE_1' ~ 5)) %>%
  ggplot(aes(y = AME, x = window_size, 
             color = specification, ymin = lower, ymax = upper)) +
  geom_pointrange(position = position_dodge(width = .4)) +
  scale_color_manual(values = c('gray10',
                                'gray30',
                                'gray50',
                                'gray70')) +
  scale_x_continuous(breaks = c(seq(1, 5, 1)), 
                     labels = c(seq(5, 1, -1))) +
  scale_y_continuous(breaks = c(seq(0, .4, .1)), 
                     labels = c(seq(0, .4, .1)), 
                     limits = c(-.01, .41)) +
  labs(y = 'Average Marginal Effect of CWE', x = 'Window Size (Years)', 
       color = '') +
  theme_pubclean() +
  theme(panel.background = element_rect(fill = 'white'),
        legend.key = element_rect(fill = 'white'),
        axis.ticks.x = element_blank()) 




# Figure 4 ----------------------------------------------------------------
#Individual minister exit probabilities by minister tenure

#Estimates use CWE_5pooled model using CWE5, full event covariates, and gwf_ssa data
#from above: cwe5[[4]][[1]]

#predicted probability of minister exit

non_cwe_period <- data.frame(minister_tenure = c(seq(1,15, .15)),
                             pers_2pl = median(gwf_ssa$pers_2pl, na.rm = T),
                             CWE_5 = 0,
                             CWE_4 = 0, 
                             CWE_3 = 0, 
                             CWE_2 = 0, 
                             CWE_1 = 0,
                             leader_tenure = median(gwf_ssa$leader_tenure, na.rm = T),
                             leader_exit = 0,
                             coup_corrected = 0,
                             election_corrected = 0,
                             cabinet_size = median(gwf_ssa$cabinet_size, na.rm = T),
                             lcgdppc = median(gwf_ssa$lcgdppc, na.rm = T),
                             ln_oilgas_rents = median(gwf_ssa$ln_oilgas_rents, na.rm = T), 
                             left_truncated = 0)

cwe_period <- data.frame(minister_tenure = c(seq(1,15, .15)),
                         pers_2pl = median(gwf_ssa$pers_2pl, na.rm = T),
                         CWE_5 = 1,
                         CWE_4 = 1, 
                         CWE_3 = 1, 
                         CWE_2 = 1, 
                         CWE_1 = 1,
                         leader_tenure = median(gwf_ssa$leader_tenure, na.rm = T),
                         leader_exit = 0,
                         coup_corrected = 0,
                         election_corrected = 0,
                         cabinet_size = median(gwf_ssa$cabinet_size, na.rm = T),
                         lcgdppc = median(gwf_ssa$lcgdppc, na.rm = T),
                         ln_oilgas_rents = median(gwf_ssa$ln_oilgas_rents, na.rm = T), 
                         left_truncated = 0)

tenure_plot <- bind_rows(non_cwe_period, cwe_period)


exit_probs <- data.frame(predict(cwe5[[4]][[1]],
                                 newdata = tenure_plot, 
                                 type = 'response', se.fit = T)) %>%
  bind_cols(., tenure_plot) %>%
  mutate(ll = fit - (1.96*se.fit),
         ul = fit + (1.96*se.fit)) 



ggplot(exit_probs[exit_probs$CWE_5==0,], aes(x = minister_tenure, y = fit)) +
  geom_line() +
  geom_ribbon(aes(ymin = ll, ymax = ul),
              fill = "grey60", alpha = .2) +
  labs(x = 'Minister Tenure', y = 'Pr(Minister Exit)') + 
  scale_x_continuous(breaks = c(1, 5, 10, 15), 
                     limits = c(1, 16)) +
  scale_y_continuous(breaks = c(seq(.1, .2, .02)), 
                     labels = c(seq(.1, .2, 0.02)), 
                     limits = c(.1, .2)) +
  theme_pubclean() 

# Figure 5 -----------------------------------------------------------
#Difference in minister exit probabilities at CWE = 1 and CWE = 0 for
#junior and senior ministers

#zelig logit function
zelig_pooled_logit <- function(xvar){
  
  zelig(as.formula(paste0('spell_end ~', xvar,'*senior', 
                          '+ ln_oilgas_rents +
                   pers_2pl +
                   leader_exit + 
                   leader_tenure +
                   I(leader_tenure^2) +
                   I(leader_tenure^3) +
                   election_corrected + 
                   coup_corrected + 
                   lcgdppc + 
                   cabinet_size + 
                   left_truncated')),
        data = gwf_ssa, 
        model = 'logit', 
        cite = F)
}


#5 year window
z.5year <- zelig_pooled_logit('CWE_5')

#4 year window
z.4year <- zelig_pooled_logit("CWE_4")

#3 year window
z.3year <- zelig_pooled_logit('CWE_3')

#2 year window
z.2year <- zelig_pooled_logit('CWE_2')

#1 year window
z.1year <- zelig_pooled_logit('CWE_1')


#First difference probability simulations

#function to extract and store simulation information
zelig_results <- function(sim_output){
  
  first_diffs <- as.numeric(unlist(sim_output$sim.out$x1$fd))
  
  data.frame(fd = mean(first_diffs), 
             std_d = sd(first_diffs), 
             ll = quantile(first_diffs, .025), 
             ul = quantile(first_diffs, .975))
  
}


#junior ministers 5 year window
x.5.noncwe_junior <- setx(z.5year, CWE_5 = 0, senior = 0)
x.5.cwe_junior <- setx(z.5year, CWE_5 = 1, senior = 0)
s.5.out_junior <- sim(z.5year, x = x.5.noncwe_junior, x1 = x.5.cwe_junior)


#junior ministers 4 year window
x.4.noncwe_junior <- setx(z.4year, CWE_4 = 0, senior = 0)
x.4.cwe_junior <- setx(z.4year, CWE_4 = 1, senior = 0)
s.4.out_junior <- sim(z.4year, x = x.4.noncwe_junior, x1 = x.4.cwe_junior)
summary(s.4.out_junior)

#junior ministers 3 year window
x.3.noncwe_junior <- setx(z.3year, CWE_3 = 0, senior = 0)
x.3.cwe_junior <- setx(z.3year, CWE_3 = 1, senior = 0)
s.3.out_junior <- sim(z.3year, x = x.3.noncwe_junior, x1 = x.3.cwe_junior)
summary(s.3.out_junior)


#junior ministers 2 year window
x.2.noncwe_junior <- setx(z.2year, CWE_2 = 0, senior = 0)
x.2.cwe_junior <- setx(z.2year, CWE_2 = 1, senior = 0)
s.2.out_junior <- sim(z.2year, x = x.2.noncwe_junior, x1 = x.2.cwe_junior)
summary(s.2.out_junior)


#junior minister 1 year window
x.1.noncwe_junior <- setx(z.1year, CWE_1 = 0, senior = 0)
x.1.cwe_junior <- setx(z.1year, CWE_1 = 1, senior = 0)
s.1.out_junior <- sim(z.1year, x = x.1.noncwe_junior, x1 = x.1.cwe_junior)
summary(s.1.out_junior)

#senior ministers 5 year window
x.5.noncwe_senior <- setx(z.5year, CWE_5 = 0, senior = 1)
x.5.cwe_senior <- setx(z.5year, CWE_5 = 1, senior = 1)
s.5.out_senior <- sim(z.5year, x = x.5.noncwe_senior, x1 = x.5.cwe_senior)
summary(s.5.out_senior)

#senior ministers 4 year window
x.4.noncwe_senior <- setx(z.4year, CWE_4 = 0, senior = 1)
x.4.cwe_senior <- setx(z.4year, CWE_4 = 1, senior = 1)
s.4.out_senior <- sim(z.4year, x = x.4.noncwe_senior, x1 = x.4.cwe_senior)
summary(s.4.out_senior)


#senior ministers 3 year window
x.3.noncwe_senior <- setx(z.3year, CWE_3 = 0, senior = 1)
x.3.cwe_senior <- setx(z.3year, CWE_3 = 1, senior = 1)
s.3.out_senior <- sim(z.3year, x = x.3.noncwe_senior, x1 = x.3.cwe_senior)
summary(s.3.out_senior)


#senior ministers 2 year window
x.2.noncwe_senior <- setx(z.2year, CWE_2 = 0, senior = 1)
x.2.cwe_senior <- setx(z.2year, CWE_2 = 1, senior = 1)
s.2.out_senior <- sim(z.2year, x = x.2.noncwe_senior, x1 = x.2.cwe_senior)
summary(s.2.out_senior)


#senior ministers 1 year window
x.1.noncwe_senior <- setx(z.1year, CWE_1 = 0, senior = 1)
x.1.cwe_senior <- setx(z.1year, CWE_1 = 1, senior = 1)
s.1.out_senior <- sim(z.1year, x = x.1.noncwe_senior, x1 = x.1.cwe_senior)


#extract all simulation results to dataframe

all_sims <- list(s.5.out_junior, s.4.out_junior, 
                 s.3.out_junior, s.2.out_junior,
                 s.1.out_junior, s.5.out_senior,
                 s.4.out_senior, s.3.out_senior, 
                 s.2.out_senior, s.1.out_senior)

all_sims_results <- lapply(all_sims, zelig_results) %>%
  bind_rows(.) %>%
  mutate(minister = c(rep("Junior", 5), rep("Senior", 5)), 
         window_size = c(rep(seq(5,1,-1), 2)), 
         window_size = factor(window_size, levels = c(5, 4, 3, 2, 1)))


#Figure 5 plot
ggplot(all_sims_results, aes(x = window_size, y = fd, ymin = ll, ymax = ul)) +
  geom_pointrange() +
  facet_wrap(~minister) +
  scale_y_continuous(breaks = c(seq(0, .6, .1)),
                     limits = c(0, .61)) +
  labs(x = 'Window Size (Years)', y = 'First Difference') +
  theme_pubclean()


# Figure 6 ----------------------------------------------------------------
# Survival Critical Portfolios

#zelig logit function
zelig_pooled_coercive_logit <- function(xvar){
  
  zelig(as.formula(paste0('spell_end ~', xvar,"*",'survival_portfolio', 
                          '+ minister_tenure +
                          I(minister_tenure^2) +
                          I(minister_tenure^3) +
                   ln_oilgas_rents +
                   pers_2pl +
                   leader_exit + 
                   leader_tenure +
                   I(leader_tenure^2) +
                   I(leader_tenure^3) +
                   election_corrected + 
                   coup_corrected + 
                   lcgdppc + 
                   cabinet_size + 
                   left_truncated')),
        data = gwf_ssa, 
        model = 'logit', 
        cite = F)
}


#function to extract and store simulation information
zelig_results <- function(sim_output){
  
  first_diffs <- as.numeric(unlist(sim_output$sim.out$x1$fd))
  
  data.frame(fd = mean(first_diffs), 
             std_d = sd(first_diffs), 
             ll = quantile(first_diffs, .025), 
             ul = quantile(first_diffs, .975))
  
}

#5 year window, high prestige
z.5year_coercive <- zelig_pooled_coercive_logit('CWE_5')
summary(z.5year_coercive)

#4 year window
z.4year_coercive <- zelig_pooled_coercive_logit("CWE_4")

#3 year window
z.3year_coercive <- zelig_pooled_coercive_logit('CWE_3')

#2 year window
z.2year_coercive <- zelig_pooled_coercive_logit('CWE_2')

#1 year window
z.1year_coercive <- zelig_pooled_coercive_logit('CWE_1')
summary(z.1year_coercive)

#High prestige/non-high prestige simulations

#non-high prestige ministers 5 year window
x.5.noncwe_nc <- setx(z.5year_coercive, CWE_5 = 0, survival_portfolio = 0)
x.5.cwe_nc <- setx(z.5year_coercive, CWE_5 = 1, survival_portfolio = 0)
s.5.out_nc <- sim(z.5year_coercive, x = x.5.noncwe_nc, x1 = x.5.cwe_nc)



#non-high prestige ministers 4 year window
x.4.noncwe_nc <- setx(z.4year_coercive, CWE_4 = 0, survival_portfolio = 0)
x.4.cwe_nc <- setx(z.4year_coercive, CWE_4 = 1, survival_portfolio = 0)
s.4.out_nc <- sim(z.4year_coercive, x = x.4.noncwe_nc, x1 = x.4.cwe_nc)
summary(s.4.out_nc)

#non-high prestige ministers 3 year window
x.3.noncwe_nc <- setx(z.3year_coercive, CWE_3 = 0, survival_portfolio = 0)
x.3.cwe_nc <- setx(z.3year_coercive, CWE_3 = 1, survival_portfolio = 0)
s.3.out_nc <- sim(z.3year_coercive, x = x.3.noncwe_nc, x1 = x.3.cwe_nc)
summary(s.3.out_nc)



#non-high prestige ministers 2 year window
x.2.noncwe_nc <- setx(z.2year_coercive, CWE_2 = 0, survival_portfolio = 0)
x.2.cwe_nc <- setx(z.2year_coercive, CWE_2 = 1, survival_portfolio = 0)
s.2.out_nc <- sim(z.2year_coercive, x = x.2.noncwe_nc, x1 = x.2.cwe_nc)
summary(s.2.out_nc)



#non-high prestige ministers 1 year window
x.1.noncwe_nc <- setx(z.1year_coercive, CWE_1 = 0, survival_portfolio = 0)
x.1.cwe_nc <- setx(z.1year_coercive, CWE_1 = 1, survival_portfolio = 0)
s.1.out_nc <- sim(z.1year_coercive, x = x.1.noncwe_nc, x1 = x.1.cwe_nc)
summary(s.1.out_nc)

#high prestige ministers 5 year window
x.5.noncwe_c <- setx(z.5year_coercive, CWE_5 = 0, survival_portfolio = 1)
x.5.cwe_c <- setx(z.5year_coercive, CWE_5 = 1, survival_portfolio = 1)
s.5.out_c <- sim(z.5year_coercive, x = x.5.noncwe_c, x1 = x.5.cwe_c)
summary(s.5.out_c)

#high prestige ministers 4 year window
x.4.noncwe_c <- setx(z.4year_coercive, CWE_4 = 0, survival_portfolio  = 1)
x.4.cwe_c <- setx(z.4year_coercive, CWE_4 = 1, survival_portfolio  = 1)
s.4.out_c <- sim(z.4year_coercive, x = x.4.noncwe_c, x1 = x.4.cwe_c)
summary(s.4.out_c)


#high prestige ministers 3 year window
x.3.noncwe_c <- setx(z.3year_coercive, CWE_3 = 0, survival_portfolio = 1)
x.3.cwe_c <- setx(z.3year_coercive, CWE_3 = 1, survival_portfolio = 1)
s.3.out_c <- sim(z.3year_coercive, x = x.3.noncwe_c, x1 = x.3.cwe_c)
summary(s.3.out_c)



#high prestige ministers 2 year window
x.2.noncwe_c <- setx(z.2year_coercive, CWE_2 = 0,survival_portfolio = 1)
x.2.cwe_c <- setx(z.2year_coercive, CWE_2 = 1, survival_portfolio = 1)
s.2.out_c <- sim(z.2year_coercive, x = x.2.noncwe_c, x1 = x.2.cwe_c)
summary(s.2.out_c)




#senior ministers 1 year window
x.1.noncwe_c <- setx(z.1year_coercive, CWE_1 = 0, survival_portfolio = 1)
x.1.cwe_c <- setx(z.1year_coercive, CWE_1 = 1, survival_portfolio = 1)
s.1.out_c <- sim(z.1year_coercive, x = x.1.noncwe_c, x1 = x.1.cwe_c)

#extract all simulation results to dataframe

all_sims_coercive <- list(s.5.out_nc, s.4.out_nc, 
                          s.3.out_nc, s.2.out_nc,
                          s.1.out_nc, s.5.out_c,
                          s.4.out_c, s.3.out_c, 
                          s.2.out_c, s.1.out_c)

all_sims_results_coercive <- lapply(all_sims_coercive, zelig_results) %>%
  bind_rows(.) %>%
  mutate(minister = c(rep("Other Portfolio", 5), rep("Survival-critical Portfolio", 5)),
         minister = factor(minister, levels = c("Other Portfolio",
                                                "Survival-critical Portfolio")),
         window_size = c(rep(seq(5,1,-1), 2)), 
         window_size = factor(window_size, levels = c(5, 4, 3, 2, 1))) 



ggplot(all_sims_results_coercive, aes(x = window_size, y = fd, ymin = ll, ymax = ul)) +
  geom_pointrange() +
  facet_wrap(~minister) +
  scale_y_continuous(breaks = c(seq(0, .6, .1)),
                     limits = c(0, .61)) +
  labs(x = 'Window Size (Years)', y = 'First Difference') +
  theme_pubclean()


# Figure 7 ----------------------------------------------------------------
# Placebo test examining all 5-year time windows

#5 year CWE placebo test function

placebo_test_5_pooled <- function(start_year){
  
  df_placebo <- gwf_ssa %>%
    mutate(time_window = ifelse(year >= start_year & year <= start_year+4, 1, 0)) 
  
  model <- glm(spell_end ~ time_window, 
               data = df_placebo, 
               family = binomial)
  
  
  robust_ses <- coeftest(model, cluster.vcov(model, df_placebo$leaderID))
  
  robust_ses <- broom::tidy(robust_ses) %>%
    filter(row_number()==2)
  
  robust_ses
}


#set start year to run through function above
start_year <- 1966:2006


cwe5_placebo_pooled <- lapply(start_year, placebo_test_5_pooled)

cwe5_placebo_pooled_results <- lapply(cwe5_placebo_pooled, `[`, c('estimate', 'std.error')) 

cwe5_placebo_pooled_results <- data.frame(matrix(unlist(cwe5_placebo_pooled_results),
                                                 nrow = length(cwe5_placebo_pooled_results), byrow = T)) %>%
  mutate(start_year = 1966:2006, 
         end_year = start_year + 4,
         ul = X1 + (1.96*X2),
         ll = X1 - (1.96*X2),
         window_used = ifelse(start_year==1988, 'Yes', 'No')) %>%
  rename(estimate = X1, 
         se = X2)

options(scipen = 999)

ggplot(cwe5_placebo_pooled_results, aes(x =start_year, y = estimate, 
                                        ymin = ll, ymax = ul, color = as.factor(window_used))) +
  geom_pointrange(size = .3) + 
  geom_hline(yintercept = 0, linetype = 'dashed') +
  scale_x_continuous(breaks = c(seq(1966, 2006, 1))) +
  scale_color_manual(values = c('gray75', 'black'), 
                     guide_legend('Chosen Window?')) +
  scale_y_continuous(limits = c(-.8, .85), 
                     breaks = c(seq(-.8, .8, .2)), 
                     labels = c(seq(-.8, .8, .2))) +
  labs(x = 'Window Start Year', y = 'Coefficient Estimate') +
  theme_pubclean() +
  theme(text = element_text(size = 12), 
        axis.text.x = element_text(angle = 75, vjust = .4),
        legend.position = 'bottom', 
        legend.key = element_blank())



# Figure 8 ----------------------------------------------------------------
# Mean number of parties in the cabinet in African autocracies

#collapse gwf_ssa data to country-year
party_df <- gwf_ssa %>%
  dplyr::select(cowcode, year, party) %>%
  filter(!is.na(party)) %>%
  group_by(cowcode, year) %>%
  summarise(parties = n_distinct(party)) %>%
  group_by(year) %>%
  summarise(parties = mean(parties))

ggplot(party_df, aes(x = year, y = parties)) +
  geom_line() +
  theme_pubclean() +
  theme(axis.text = element_text(size = 12), 
        axis.title = element_text(size = 12), 
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12), 
        legend.position = 'bottom', 
        legend.background = element_blank(),
        legend.key = element_blank()) +
  labs(x = 'Year', y = '# of Parties in Cabinet')  



# Figure 9 ----------------------------------------------------------------
# Probability of Senior Minister Survival during Cold War,
#Cold War End, and post-Cold War

gwf_ssa <- gwf_ssa %>%
  mutate(PCW = ifelse(year > 1992, 1, 0))

senior_cwe5_pcw <- glm(spell_end ~ CWE_5*senior +
                         PCW*senior +
                         ln_oilgas_rents +
                         pers_2pl +
                         leader_exit +
                         leader_tenure +
                         I(leader_tenure^2) +
                         I(leader_tenure^3) +
                         election_corrected +
                         coup_corrected +
                         lcgdppc +
                         cabinet_size +
                         left_truncated,
                       data = gwf_ssa, 
                       family = binomial)

senior_cwe5_pcw_se <- sqrt(diag(cluster.vcov(senior_cwe5_pcw,
                                             gwf_ssa$leaderID)))


pred_df <- data.frame(CWE_5 = c(0, 1, 0, 0, 1, 0),
                      PCW = c(0, 0, 1, 0, 0, 1),
                      senior = c(0, 0 , 0, 1, 1, 1),
                      minister_tenure = c(rep(median(gwf_ssa$minister_tenure, na.rm = T), 6)), 
                      ln_oilgas_rents = c(rep(median(gwf_ssa$ln_oilgas_rents, na.rm = T), 6)),
                      pers_2pl = c(rep(median(gwf_ssa$pers_2pl, na.rm = T), 6)), 
                      leader_exit = c(rep(0, 6)), 
                      leader_tenure = c(rep(median(gwf_ssa$leader_tenure, na.rm = T), 6)), 
                      election_corrected = c(rep(0, 6)), 
                      coup_corrected = c(rep(0, 6)), 
                      lcgdppc = c(rep(median(gwf_ssa$lcgdppc, na.rm = T), 6)),
                      cabinet_size = c(rep(median(gwf_ssa$cabinet_size, na.rm = T), 6)), 
                      left_truncated = c(rep(0, 6)))

senior_cwe5_pcw_pred <- data.frame(predict(senior_cwe5_pcw, newdata = pred_df, 
                                           se.fit = T, 
                                           type = 'response')) %>%
  mutate( fit = 1 - fit, 
          ll = fit - (1.96*se.fit), 
          ul = fit + (1.96*se.fit), 
          period = c('Cold War', 'CWE', 'Post-Cold War', 
                     'Cold War', 'CWE', 'Post-Cold War'), 
          minister = c(rep('Junior', 3), 
                       rep('Senior', 3))) %>%
  filter(minister == 'Senior')




ggplot(senior_cwe5_pcw_pred, aes(x = period, year, y = fit, 
                                 ymin = ll, ymax = ul)) +
  geom_pointrange() +
  labs(x = '', y = 'Pr(Minister Survival)') +
  scale_y_continuous(breaks = c(.75, .80, .85, .90)) +
  theme_pubclean() +
  theme(text = element_text(size = 12))


# Supplemental Material Table A1 ------------------------------------------

#function for displaying country-years in sample
cy_table <- function(x){
  
  x %>%
    group_by(country_name) %>%
    #identify breaks in authoritarianism
    mutate(year_l = dplyr::lag(year, n = 1L), 
           year_l = ifelse(is.na(year_l), year - 1, year_l),
           gap = ifelse(year - year_l > 1, 1, 0), 
           spell = cumsum(gap)) %>%
    group_by(country_name, spell) %>%
    summarise(country_name = first(country_name), 
              start_year = min(year),
              end_year = max(year)) %>%
    dplyr::select(-spell)
}

#summary statistics function
summary_table <- function(x){
  x %>%
    dplyr::select(spell_end, 
                  CWE_5, 
                  minister_tenure,
                  senior,
                  high_prestige,
                  leader_tenure, 
                  leader_exit, 
                  pers_2pl, 
                  election_corrected,
                  coup_corrected,
                  lcgdppc,
                  ln_oilgas_rents, 
                  ln_total_aid_AD,
                  v2x_polyarchy, 
                  cabinet_size,
                  left_truncated)
}


# Main GWF SSA Sample Summary Statistics 

gwf_ssa_table <- cy_table(gwf_ssa)


# GWF All Africa Sample Summary Statistics

gwf_all <- df %>%
  filter(!is.na(gwf_regimetype))

gwf_all_table <- cy_table(gwf_all)


# Lexical Index of Electoral Democracy sub-Saharan Sample Summary Statistics

lied_ssa <- df %>%
  filter(lied_autoc==1) %>%
  filter(!(country_name %in% nafrica_sa))

lied_ssa_table <- cy_table(lied_ssa) 

# Lexical Index of Electoral Democracy All Africa Sample Summary Statistics

lied_all <- df %>%
  filter(lied_autoc==1) 

lied_all_table <- cy_table(lied_all)

# VDem ROW sub-Saharan Sample Summary Statistics 

row_ssa <- df %>%
  filter(vdem_autoc==1) %>%
  filter(!(country_name %in% nafrica_sa))

row_ssa_table <- cy_table(row_ssa)

# VDem ROW All Africa Sample Summary Statistics 

row_all <- df %>%
  filter(vdem_autoc==1) 

row_all_table <- cy_table(row_all)

# Balanced sub-Saharan Sample Summary Statistics 

#find continuous cases of autocracy using gwf_ssa sample 
#cwe start = 1988, cwe end = 1992

ssa_balance_remove <- c('Morocco','Algeria','Tunisia','Libya','Egypt','South Africa',
                        'Namibia', 'Somalia')

gwf_ssa_cw_regimes <- df %>%
  filter(!(country_name %in% ssa_balance_remove)) %>%
  mutate(cold_war = ifelse(year < 1988, 1, 0), 
         gwf_dictatorship = ifelse(!is.na(gwf_regimetype), 1, 0)) %>%
  filter(cold_war==1) %>%
  group_by(cowcode) %>%
  mutate(dictatorship_percent = mean(gwf_dictatorship)) 

gwf_always <- gwf_ssa_cw_regimes %>%
  filter(dictatorship_percent==1) %>%
  group_by(cowcode) %>%
  summarize(country_name = first(country_name)) %>%
  pull(cowcode)

gwf_ssa_balanced <- df %>%
  #select countries that had uninterrupted autocratic spells (not necessarily same regime)
  #going into the cwe period
  filter(cowcode %in% gwf_always & cowcode != 553| #malawi excluded because of missing 1983 data
           #benin post 1972
           cowcode==434 & year > 1972| 
           #sierra leone post 1966
           cowcode==451 & year > 1966| 
           #ghana post 1981
           cowcode==452 & year > 1981| 
           #nigeria post 1983
           cowcode==475 & year > 1983| 
           #chad post 1981
           cowcode==483 & year > 1981| 
           #uganda post 1980
           cowcode==500 & year > 1980| 
           #angola post 1975
           cowcode==540 & year > 1975| 
           #zambia post 1966
           cowcode==551 & year > 1966| 
           #lesotho after 1969
           cowcode==570 & year > 1969) %>%
  filter(year < 2011)



# Country-Years and Minister-Years in All Samples 

#list of all samples
samples <- list(gwf_ssa, gwf_all, 
                lied_ssa, lied_all,
                row_ssa, row_all)


cy_snapshot <- function(x){
  
  country_years <- x %>%
    dplyr::select(cowcode, year) %>%
    group_by(cowcode, year) %>%
    summarise_all(first) %>%
    ungroup() %>%
    summarise(Countries = n_distinct(cowcode), 
              `Country-years` = n())
  
  minister_years <- x %>%
    dplyr::select(ministerID, year) %>%
    group_by(ministerID, year) %>%
    summarise_all(first) %>%
    ungroup() %>%
    summarise(Ministers = n_distinct(ministerID),
              `Minister-years` = n())
  
  bind_cols(country_years, minister_years)
}

samples_snapshot <- lapply(samples, cy_snapshot) %>%
  bind_rows(.) %>%
  mutate(Sample = c('GWF sub-Saharan (Main)', 'GWF All', 
                    'LIED sub-Saharan', 'LIED All', 
                    'VDem ROW sub-Saharan', 'VDem ROW All'))%>%
  dplyr::select(Sample, Countries, `Country-years`, 
                Ministers, `Minister-years`)

#Table A1
stargazer(samples_snapshot, 
          summary = F, 
          rownames = F,
          table.placement = '!hbt',
          title = 'Summary of Countries and Ministers across Sample Variations',
          out = 'samples_snapshot.tex')





# Supplemental Material Table A2 ------------------------------------------

gwf_ssa_summary <- summary_table(gwf_ssa)

stargazer(gwf_ssa_summary, 
          covariate.labels = c('Minister Exit', 
                               'CWE 5-year', 
                               'Minister Tenure',
                               'Senior Minister',
                               'High Prestige Portfolio',
                               'Leader Tenure', 
                               'Leader Exit',
                               'Personalism',
                               'Election', 
                               'Failed Coup',
                               'Ln(GDP per Capita)',
                               'Cabinet Size',
                               'Ln(Oil/Gas Rents)',
                               'Ln(Total Aid)',
                               'VDem Polyarchy',
                               'Left Truncated'),
          table.placement = '!hbt',
          summary.stat = c('n', 'mean', 'sd', 'min', 'p25', 'p75', 'max'),
          digits = 2,
          title = 'GWF Autocracy, sub-Saharan Sample Summary Statistics', 
          out = 'gwf_ssa_stats.tex')

# Supplemental Material Table A3 ------------------------------------------


gwf_ssa_table <- cy_table(gwf_ssa)

print(xtable(gwf_ssa_table, 
             caption = 'GWF Autocracy, sub-Saharan Sample Country-Years'),
      include.rownames = FALSE,
      type = 'latex',
      tabular.environment = 'longtable',
      floating = FALSE,
      caption.placement = 'top', 
      file = 'gwf_ssa_sample.tex')


# Supplemental Material Table A4 ------------------------------------------

gwf_all_summary <- summary_table(gwf_all)

stargazer(gwf_all_summary, 
          covariate.labels = c('Minister Exit', 
                               'CWE 5-year', 
                               'Minister Tenure',
                               'Senior Minister',
                               'High Prestige Portfolio',
                               'Leader Tenure', 
                               'Leader Exit',
                               'Personalism',
                               'Election', 
                               'Failed Coup',
                               'Ln(GDP per Capita)',
                               'Cabinet Size',
                               'Ln(Oil/Gas Rents)',
                               'Ln(Total Aid)',
                               'VDem Polyarchy',
                               'Left Truncated'),
          table.placement = '!hbt',
          summary.stat = c('n', 'mean', 'sd', 'min', 'p25', 'p75', 'max'),
          digits = 2,
          title = 'GWF Autocracy, All Africa Sample Summary Statistics', 
          out = 'gwf_all_stats.tex')


# Supplemental Material Table A5 ------------------------------------------

gwf_all_table <- cy_table(gwf_all)

print(xtable(gwf_all_table, 
             caption = 'GWF Autocracy, All Africa Sample Country-Years'),
      include.rownames = FALSE,
      type = 'latex',
      tabular.environment = 'longtable',
      floating = FALSE,
      caption.placement = 'top', 
      file = 'gwf_all_sample.tex')


# Supplemental Material A14 -----------------------------------------------


gwf_ssa_balanced_table <- cy_table(gwf_ssa_balanced)

print(xtable(gwf_ssa_balanced_table, 
             caption = 'GWF Autocracy, sub-Saharan Balanced Sample'),
      include.rownames = FALSE,
      type = 'latex',
      tabular.environment = 'longtable',
      floating = FALSE,
      caption.placement = 'top', 
      file = 'gwf_ssa_balanced_sample.tex')


# Supplemental Material Table A15 -----------------------------------------

gwf_ssa_balanced_summary <- summary_table(gwf_ssa_balanced)

stargazer(gwf_ssa_balanced_summary, 
          covariate.labels = c('Minister Exit', 
                               'CWE 5-year', 
                               'Minister Tenure',
                               'Senior Minister',
                               'High Prestige Portfolio',
                               'Leader Tenure', 
                               'Leader Exit',
                               'Personalism',
                               'Election', 
                               'Failed Coup',
                               'Ln(GDP per Capita)',
                               'Cabinet Size',
                               'Ln(Oil/Gas Rents)',
                               'Ln(Total Aid)',
                               'VDem Polyarchy',
                               'Left Truncated'),
          table.placement = '!hbt',
          digits = 2,
          title = 'GWF Autocracy, sub-Saharan Balanced Sample Summary Statistics', 
          out = 'gwf_ssa_balanced_stats.tex')

# Supplemental Material Table A16 -----------------------------------------

#find cwe surviving leaders
survivors <- gwf_ssa %>%
  dplyr::select(leader_name, leaderID, country_name, year, CWE_5) %>%
  group_by(leaderID, year) %>%
  summarise_all(first) %>%
  group_by(leaderID) %>%
  mutate(leader_start = first(year), 
         leader_end = last(year)) %>%
  #select only leaders with at least three years of data before and
  #after CWE window. Any less will probably not satisfy reviewers.
  filter(leader_start < 1985,
         leader_end > 1995) %>%
  summarise_all(first)  


survivors_table <- survivors %>%
  dplyr::select(country_name, leader_name, leader_start, leader_end) %>%
  rename(Country = country_name, 
         Leader = leader_name, 
         Start = leader_start, 
         End = leader_end) %>%
  arrange(Country)

stargazer(survivors_table, 
          summary = F,
          rownames = F,
          title = 'CWE Surviving Leaders, GWF SSA Sample',
          notes = 'Note: Leader start and end years reflect cabinet observations in our dataset.',
          out = 'cwe_survivors.tex')


# Supplemental Material Table A17 -----------------------------------------

#function to apply pooled logit models across data samples

samples_pooled_logit <- function(xvar, controls){
  
  lapply(samples, function(y){
    
    model <- glm(as.formula(paste('spell_end ~', 
                                  xvar, controls)), 
                 data = y, 
                 family = binomial)
    
    
    SEs <- sqrt(diag(cluster.vcov(model, y$leaderID)))
    
    return(list(model, SEs))
    
  })
}


#Table A17, no controls
A17 <- samples_pooled_logit('CWE_5', none)


#no controls table of results
stargazer(A17[[1]][[1]], A17[[2]][[1]], 
          A17[[3]][[1]], A17[[4]][[1]], 
          A17[[5]][[1]], A17[[6]][[1]], 
          se = list(A17[[1]][[2]], A17[[2]][[2]], 
                    A17[[3]][[2]], A17[[4]][[2]], 
                    A17[[5]][[2]], A17[[6]][[2]]), 
          table.placement = '!hbt',
          dep.var.labels = 'Minister Exit',
          column.labels = ' &  & &  & VDem & VDem \\\\
& GWF & GWF & LIED & LIED & ROW & ROW \\\\
&    SSA     &  All     &    SSA    &   All  &     SSA     & All\\\\',
          covariate.labels = 'CWE 5-year',
          omit.stat = c('ll', 'aic'),
          digits = 2,
          no.space = T,
          title = 'Pooled Logit of Minister Survival with 5-year CWE Window, No Controls', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          label = 'tab:cwe5_pooled_none',
          out = 'cwe5_pooled_none.tex')



# Supplemental Material Table A18 -----------------------------------------

A18 <- samples_pooled_logit('CWE_5', baseline)

stargazer(A18[[1]][[1]], A18[[2]][[1]], 
          A18[[3]][[1]], A18[[4]][[1]], 
          A18[[5]][[1]], A18[[6]][[1]], 
          se = list(A18[[1]][[2]], A18[[2]][[2]], 
                    A18[[3]][[2]], A18[[4]][[2]], 
                    A18[[5]][[2]], A18[[6]][[2]]), 
          table.placement = '!hbt',
          dep.var.labels = 'Minister Exit',
          column.labels = '  &  & &  & VDem & VDem \\\\
& GWF & GWF & LIED & LIED & ROW & ROW \\\\
&    SSA     &  All     &    SSA    &   All  &     SSA     & All\\\\',
          label = 'tab:cwe5_pooled_baseline',
          covariate.labels = c('CWE 5-year', 'Ln(Oil/Gas Rents)',
                               'Personalism', 'Ln(GDP per Capita)', 
                               'Cabinet Size', 'Left Truncated'),
          omit.stat = c('ll', 'aic'),
          digits = 2,
          no.space = TRUE,
          title = 'Pooled Logit of Minister Survival with 5-year CWE Window, Baseline Controls', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          out = 'cwe5_pooled_baseline.tex')


# Supplemental Material Table A19 -----------------------------------------


A19 <- samples_pooled_logit('CWE_5', tenure)

stargazer(A19[[1]][[1]], A19[[2]][[1]], 
          A19[[3]][[1]], A19[[4]][[1]], 
          A19[[5]][[1]], A19[[6]][[1]], 
          se = list(A19[[1]][[2]], A19[[2]][[2]], 
                    A19[[3]][[2]], A19[[4]][[2]], 
                    A19[[5]][[2]], A19[[6]][[2]]), 
          table.placement = '!hbt',
          
          dep.var.labels = 'Minister Exit',
          column.labels = '  &  & &  & VDem & VDem \\\\
& GWF & GWF & LIED & LIED & ROW & ROW \\\\
&    SSA     &  All     &    SSA    &   All  &     SSA     & All\\\\',
          label = 'tab:cwe5_pooled_tenure',
          omit.stat = c('ll', 'aic'),
          digits = 2,
          covariate.labels = c('CWE 5-year', 'Minister Tenure',
                               'Minister Tenure$^{2}$', 
                               'Minister Tenure$^{3}$',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 'Left Truncated'),
          no.space = T,
          title = 'Pooled Logit of Minister Survival with 5-year CWE Window, Tenure Controls', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          out = 'cwe5_pooled_tenure.tex')



# Supplemental Material Table A20 -----------------------------------------

A20 <- samples_pooled_logit('CWE_5', events)

stargazer(A20[[1]][[1]], A20[[2]][[1]], 
          A20[[3]][[1]], A20[[4]][[1]], 
          A20[[5]][[1]], A20[[6]][[1]], 
          se = list(A20[[1]][[2]], A20[[2]][[2]], 
                    A20[[3]][[2]], A20[[4]][[2]], 
                    A20[[5]][[2]], A20[[6]][[2]]), 
          table.placement = '!hbt',
          
          dep.var.labels = 'Minister Exit',
          column.labels = '  &  & &  & VDem & VDem \\\\
& GWF & GWF & LIED & LIED & ROW & ROW \\\\
&    SSA     &  All     &    SSA    &   All  &     SSA     & All\\\\',
          label = 'tab:cwe5_pooled_events',
          omit.stat = c('ll', 'aic'),
          digits = 2,
          covariate.labels = c('CWE 5-year', 'Minister Tenure',
                               'Minister Tenure$^{2}$', 
                               'Minister Tenure$^{3}$',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 'Leader Exit',
                               'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Election',
                               'Failed Coup',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 'Left Truncated'),
          no.space = T,
          title = 'Pooled Logit of Minister Survival with 5-year CWE Window, Events Controls', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          out = 'cwe5_pooled_events.tex')



# Supplemental Material Table A27 -----------------------------------------


pooled_balanced_logit <- function(xvar, controls){
  
  model <- glm(as.formula(paste('spell_end ~', 
                                xvar, controls)), 
               data = gwf_ssa_balanced, 
               family = binomial)
  
  SEs <- sqrt(diag(cluster.vcov(model, gwf_ssa_balanced$leaderID)))
  
  return(list(model, SEs))
  
}

#Table A27 models
#no controls
A27_none <- pooled_balanced_logit('CWE_5', none)

#baseline controls
A27_base <- pooled_balanced_logit('CWE_5', baseline)

#tenure controls
A27_tenure <- pooled_balanced_logit('CWE_5', tenure)

#events controls
A27_events <- pooled_balanced_logit('CWE_5', events)

#Table A27
stargazer(A27_none[[1]],
          A27_base[[1]],
          A27_tenure[[1]],
          A27_events[[1]],
          se = list(A27_none[[2]],
                    A27_base[[2]],
                    A27_tenure[[2]],
                    A27_events[[2]]),
          dep.var.labels.include = F,
          dep.var.caption = '',
          column.labels = c('No Controls',
                            'Baseline',
                            'Tenure',
                            'Events'),
          label = 'tab:cwe5_pooled_balanced',
          table.placement = '!hbt',
          covariate.labels = c('CWE 5-year', 'Minister Tenure',
                               'Minister Tenure$^{2}$', 
                               'Minister Tenure$^{3}$',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 'Leader Exit',
                               'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Election',
                               'Failed Coup',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 'Left Truncated'),
          omit.stat =  c('ll', 'aic', 'rsq', 'max.rsq', 'logrank', 
                         'lr', 'wald'),
          digits = 2,
          no.space = T,
          title = 'Pooled Models of Minister Survival with 5-year CWE Window, Balanced Sample', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = 'cwe5_pooled_balanced.tex')




# Supplemental Material Table A29 -----------------------------------------


#save surviving leader ids to filter gwf_ssa data
survivor_ids <- survivors$leaderID

gwf_ssa_survivors <- gwf_ssa %>%
  filter(leaderID %in% survivor_ids)


pooled_logit_survivors <- function(controls, data_sample){
  
  model <- glm(as.formula(paste('spell_end ~ CWE_5', controls)), 
               data = data_sample, 
               family = binomial)
  
  
  SEs <- sqrt(diag(cluster.vcov(model, data_sample$leaderID)))
  
  clustered_vcov <- cluster.vcov(model, data_sample$leaderID)
  
  margins_df <- data.frame(summary(margins(model, vcov = clustered_vcov)))
  
  return(list(model, SEs, margins_df))
  
}

#Table A29 models
#pooled model with no controls
A29_none <- pooled_logit_survivors(none, gwf_ssa_survivors)

#pooled model with baseline controls
A29_baseline <- pooled_logit_survivors(baseline, gwf_ssa_survivors)

#pooled model with tenure controls
A29_tenure <- pooled_logit_survivors(tenure, gwf_ssa_survivors)

#pooled model with event controls
A29_events <- pooled_logit_survivors(events, gwf_ssa_survivors)

#Table A29
#results for survivor pooled models
stargazer(A29_none[[1]], A29_baseline[[1]], 
          A29_tenure[[1]], A29_events[[1]], 
          se = list(A29_none[[2]], 
                    A29_baseline[[2]], 
                    A29_tenure[[2]], 
                    A29_events[[2]]), 
          table.placement = '!hbt',
          dep.var.labels = 'Minister Exit',
          column.labels = c('No Controls', 'Baseline', 
                            'Tenure', 'Events'),
          covariate.labels = c('CWE 5-year', 'Minister Tenure',
                               'Minister Tenure$^{2}$', 
                               'Minister Tenure$^{3}$',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 'Leader Exit',
                               'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Election',
                               'Failed Coup',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 'Left Truncated'),
          no.space = T,
          omit.stat = c('ll', 'aic', 'rsq', 'max.rsq', 'logrank', 
                        'lr', 'wald'),
          title = 'Pooled Logit of Minister Survival, GWF SSA Sample Limited to CWE Surviving Leaders', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          out = 'survivors_pooled.tex')


# Supplemental Material Table A31 -----------------------------------------

#Pooled logit models with senior dummy by CWE window interaction
senior_pooled_logit <- function(xvar){
  
  senior_df <- gwf_ssa %>%
    mutate(cwe_option = !!sym(xvar))
  
  glm(spell_end ~ cwe_option*senior + 
        ln_oilgas_rents +
        pers_2pl +
        leader_exit + 
        leader_tenure +
        I(leader_tenure^2) +
        I(leader_tenure^3) +
        election_corrected + 
        coup_corrected + 
        lcgdppc + 
        cabinet_size + 
        left_truncated,
      data = senior_df, 
      family = binomial)
}

#Table A31 models
A31 <- lapply(cwe_options, senior_pooled_logit)


#Table A31
stargazer(A31[[1]], A31[[2]],
          A31[[3]], A31[[4]],
          A31[[5]], 
          table.placement = '!hbt',
          dep.var.labels = 'Minister Exit',
          column.labels  = c('CWE 5-year', 'CWE 4-year', 
                             'CWE 3-year', 'CWE 2-year',
                             'CWE 1-year'), 
          order = c(1, 2, 14, 3, 4, 5, 6, 7, 
                    8, 9, 10, 11, 12, 13),
          digits = 2,
          covariate.labels = c('CWE Window', 
                               'Senior Minister',
                               'CWE X Senior',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 
                               'Leader Exit',
                               'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Election',
                               'Failed Coup',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 
                               'Left Truncated'),
          no.space = T,
          omit.stat = c('ll', 'aic', 'rsq', 'max.rsq', 'logrank', 
                        'lr', 'wald'),
          title = 'Pooled Logit Models with Senior Minister by CWE Interaction from Figure 5', 
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          out = 'cwe_senior_pooled.tex')




# Supplemental Material Figure A1 -----------------------------------------


gwf_ssa_senior <- gwf_ssa %>%
  filter(senior==1)

#senior model formulas

tenure_senior <- ' + ln_oilgas_rents + 
                    pers_2pl +
                    leader_tenure +
                    I(leader_tenure^2) +
                    I(leader_tenure^3) +
                    lcgdppc + 
                    cabinet_size + 
                    left_truncated'


events_senior <- '+ ln_oilgas_rents + 
                  pers_2pl +
                  leader_exit + 
                  leader_tenure +
                  I(leader_tenure^2) +
                  I(leader_tenure^3) +
                  election_corrected + 
                  coup_corrected + 
                  lcgdppc + 
                  cabinet_size + 
                  left_truncated'


aid_dem_senior <- '+ ln_oilgas_rents +
                  ln_total_aid_AD +
                  v2x_polyarchy +
                  pers_2pl +
                  leader_exit + 
                  leader_tenure +
                  I(leader_tenure^2) +
                  I(leader_tenure^3) +
                  election_corrected + 
                  coup_corrected + 
                  lcgdppc + 
                  cabinet_size + 
                  left_truncated'

model_specs_senior <- list(none, baseline, tenure_senior, 
                           events_senior, aid_dem_senior)

Figure_A1<- lapply(cwe_options, function(x){
  
  lapply(model_specs_senior, function(controls){
    model <-  glm(as.formula(paste('spell_end ~', 
                                   x, controls)), 
                  data = gwf_ssa_senior, 
                  family = binomial)
    
    clustered_vcov <- cluster.vcov(model, gwf_ssa_senior$leaderID)
    
    #marginal effects 
    data.frame(summary(margins(model, vcov = clustered_vcov))) %>%
      mutate(model = 'Pooled', 
             sample = 'Senior Ministers') %>%
      filter(grepl('CWE', factor)==T)
  })
})


senior_margins_df <- bind_rows(Figure_A1) %>%
  mutate(specification = rep(c('No Controls', 
                               '+ Baseline',
                               '+ Tenure',
                               '+ Events',
                               '+ Extras'), 5), 
         specification = factor(specification, 
                                levels = c('No Controls', 
                                           '+ Baseline',
                                           '+ Tenure',
                                           '+ Events',
                                           '+ Extras'))) %>%
  filter(specification != '+ Extras') %>%
  mutate(window_size = case_when(factor=='CWE_5' ~ 1, 
                                 factor=='CWE_4' ~ 2, 
                                 factor=='CWE_3' ~ 3, 
                                 factor=='CWE_2' ~ 4, 
                                 factor=='CWE_1' ~ 5))


ggplot(senior_margins_df, aes(y = AME, x = window_size, 
                              color = specification, ymin = lower, ymax = upper)) +
  geom_pointrange(position = position_dodge(width = .4)) +
  scale_color_manual(values = c('gray10',
                                'gray30',
                                'gray50',
                                'gray70')) +
  scale_x_continuous(breaks = c(seq(1, 5, 1)), 
                     labels = c(seq(5, 1, -1))) +
  scale_y_continuous(breaks = c(seq(0, .4, .1)), 
                     labels = c(seq(0, .4, .1)), 
                     limits = c(-.01, .45)) +
  labs(y = 'Average Marginal Effect of CWE', x = 'Window Size (Years)', 
       color = '') +
  theme_pubclean() +
  theme(panel.background = element_rect(fill = 'white'),
        legend.key = element_rect(fill = 'white'),
        axis.ticks.x = element_blank()) 

ggsave('Fig2_Senior_AMEs.pdf', 
       width = 6, height = 4)


# Supplemental Material Table A42 -----------------------------------------


#coercive portfolio logit function
coercive_pooled_logit <- function(xvar){
  
  coercive_df <- gwf_ssa %>%
    mutate(cwe_option = !!sym(xvar))
  
  glm(spell_end ~ cwe_option*survival_portfolio + 
        ln_oilgas_rents +
        pers_2pl +
        leader_exit + 
        leader_tenure +
        I(leader_tenure^2) +
        I(leader_tenure^3) +
        election_corrected + 
        coup_corrected + 
        lcgdppc + 
        cabinet_size + 
        left_truncated,
      data = coercive_df, 
      family = binomial)
}

# Table A42 models
A42 <- lapply(cwe_options, coercive_pooled_logit)

#Table A42
stargazer(A42[[1]], A42[[2]],
          A42[[3]], A42[[4]],
          A42[[5]], 
          table.placement = '!hbt',
          dep.var.labels = 'Minister Exit',
          column.labels  = c('CWE 5-year', 'CWE 4-year', 
                             'CWE 3-year', 'CWE 2-year',
                             'CWE 1-year'), 
          order = c(1, 2, 14, 3, 4, 5, 6, 7, 
                    8, 9, 10, 11, 12, 13),
          digits = 2,
          covariate.labels = c('CWE Window', 
                               'Coercive Portfolio',
                               'CWE X Coercive Portfolio',
                               'Ln(Oil/Gas Rents)',
                               'Personalism', 
                               'Leader Exit',
                               'Leader Tenure',
                               'Leader Tenure$^{2}$', 
                               'Leader Tenure$^{3}$',
                               'Election',
                               'Failed Coup',
                               'Ln(GDP per Capita)', 
                               'Cabinet Size', 
                               'Left Truncated'),
          no.space = T,
          omit.stat = c('ll', 'aic', 'rsq', 'max.rsq', 'logrank', 
                        'lr', 'wald'),
          title = 'Pooled Logit Models with Coercive Portfolio by CWE Interaction from Figure 6', 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          label = 'tab:cwe_coercive_pooled',
          out = 'cwe_coercive_pooled.tex')


# Supplemental Material Figure A2 -----------------------------------------

gwf_ssa_coercive <- gwf_ssa %>%
  filter(survival_portfolio==1)

Figure_A2 <- lapply(cwe_options, function(x){
  
  lapply(model_specs, function(controls){
    model <-  glm(as.formula(paste('spell_end ~', 
                                   x, controls)), 
                  data = gwf_ssa_coercive, 
                  family = binomial)
    
    clustered_vcov <- cluster.vcov(model, gwf_ssa_coercive$leaderID)
    
    #marginal effects 
    data.frame(summary(margins(model, vcov = clustered_vcov))) %>%
      mutate(model = 'Pooled', 
             sample = 'Coercive Portfolio') %>%
      filter(grepl('CWE', factor)==T)
  })
})


high_margins_df <- bind_rows(Figure_A2) %>%
  mutate(specification = rep(c('No Controls', 
                               '+ Baseline',
                               '+ Tenure',
                               '+ Events'), 5), 
         specification = factor(specification, 
                                levels = c('No Controls', 
                                           '+ Baseline',
                                           '+ Tenure',
                                           '+ Events'))) %>%
  filter(specification != '+ Extras') %>%
  mutate(window_size = case_when(factor=='CWE_5' ~ 1, 
                                 factor=='CWE_4' ~ 2, 
                                 factor=='CWE_3' ~ 3, 
                                 factor=='CWE_2' ~ 4, 
                                 factor=='CWE_1' ~ 5))


ggplot(high_margins_df, aes(y = AME, x = window_size, 
                            color = specification, ymin = lower, ymax = upper)) +
  geom_pointrange(position = position_dodge(width = .4)) +
  scale_color_manual(values = c('gray10',
                                'gray30',
                                'gray50',
                                'gray70')) +
  scale_x_continuous(breaks = c(seq(1, 5, 1)), 
                     labels = c(seq(5, 1, -1))) +
  scale_y_continuous(breaks = c(seq(0, .4, .1)), 
                     labels = c(seq(0, .4, .1)), 
                     limits = c(-.01, .41)) +
  labs(y = 'Average Marginal Effect of CWE', x = 'Window Size (Years)', 
       color = '') +
  theme_pubclean() +
  theme(panel.background = element_rect(fill = 'white'),
        legend.key = element_rect(fill = 'white'),
        axis.ticks.x = element_blank()) 

ggsave('Fig2_Coercive_Portfolio_AMEs.pdf', 
       width = 6, height = 4)



# Supplemental Material Figure A6 -----------------------------------------

placebo_test_5_fe <- function(start_year){
  
  df_placebo <- gwf_ssa %>%
    mutate(time_window = ifelse(year >= start_year & year <= start_year+4, 1, 0)) 
  
  model <- clogit(spell_end ~ time_window + strata(leaderID), 
                  data = df_placebo, 
                  method = 'efron', 
                  cluster = leaderID)
  
  results <- broom::tidy(model) %>%
    dplyr::select(estimate, robust.se)
  
  results 
}


cwe5_placebo_fe <- lapply(start_year, placebo_test_5_fe)

cwe5_placebo_fe_results <- bind_rows(cwe5_placebo_fe) %>%
  mutate(ll = estimate - (1.96*robust.se), 
         ul = estimate + (1.96*robust.se), 
         start_year: 1966:2006,
         window_used = ifelse(start_year==1988, 'Yes', 'No'))


ggplot(cwe5_placebo_fe_results, aes(x = start_year, y = estimate, 
                                    ymin = ll, ymax = ul,
                                    color = as.factor(window_used))) +
  geom_pointrange(size = .3) + 
  geom_hline(yintercept = 0, linetype = 'dashed') +
  scale_x_continuous(breaks = c(seq(1966, 2006, 1))) +
  scale_color_manual(values = c('gray75', 'black'), 
                     guide_legend('Chosen Window?')) +
  scale_y_continuous(limits = c(-.8, .85), 
                     breaks = c(seq(-.8, .8, .2)), 
                     labels = c(seq(-.8, .8, .2))) +
  labs(x = 'Window Start Year', y = 'Coefficient Estimate') +
  theme_pubclean() +
  theme(text = element_text(size = 12), 
        axis.text.x = element_text(angle = 75, vjust = .4),
        legend.position = 'bottom', 
        legend.key = element_blank())

ggsave('Fig5_CWE5_placebo_fe.pdf', 
       width = 9, height = 6)



# Supplemental Material Figure A7 -----------------------------------------


placebo_test_3_pooled <- function(start_year){
  
  df_placebo <- gwf_ssa %>%
    mutate(time_window = ifelse(year >= start_year & year <= start_year+2, 1, 0)) 
  
  model <- glm(spell_end ~ time_window, 
               data = df_placebo, 
               family = binomial)
  
  
  robust_ses <- coeftest(model, cluster.vcov(model, df_placebo$leaderID))
  
  robust_ses <- broom::tidy(robust_ses) %>%
    filter(row_number()==2)
  
  robust_ses
}

start_year3 <- 1966:2008


cwe3_placebo_pooled <- lapply(start_year3, placebo_test_3_pooled)

cwe3_placebo_pooled_results <- lapply(cwe3_placebo_pooled, `[`, c('estimate', 'std.error')) 

cwe3_placebo_pooled_results <- data.frame(matrix(unlist(cwe3_placebo_pooled_results),
                                                 nrow = length(cwe3_placebo_pooled_results), byrow = T)) %>%
  mutate(start_year = 1966:2008, 
         end_year = start_year + 4,
         ul = X1 + (1.96*X2),
         ll = X1 - (1.96*X2)) %>%
  rename(estimate = X1, 
         se = X2)

ggplot(cwe3_placebo_pooled_results, aes(x =start_year, y = estimate, 
                                        ymin = ll, ymax = ul)) +
  geom_pointrange(size = .3) + 
  geom_hline(yintercept = 0, linetype = 'dashed') +
  scale_x_continuous(breaks = c(seq(1966, 2008, 1))) +
  scale_y_continuous(limits = c(-.8, 1), 
                     breaks = c(seq(-.8, 1, .2)), 
                     labels = c(seq(-.8, 1, .2))) +
  labs(x = 'Window Start Year', y = 'Coefficient Estimate') +
  theme_pubclean() +
  theme(text = element_text(size = 12), 
        axis.text.x = element_text(angle = 75, vjust = .4),
        legend.position = 'bottom', 
        legend.key = element_blank())

ggsave('Fig5_CWE3_placebo.pdf', 
       width = 9, height = 6)


# Supplemental Material Figure A8 -----------------------------------------

placebo_test_3_fe <- function(start_year){
  
  df_placebo <- gwf_ssa %>%
    mutate(time_window = ifelse(year >= start_year & year <= start_year+2, 1, 0)) 
  
  model <- clogit(spell_end ~ time_window + strata(leaderID), 
                  data = df_placebo, 
                  method = 'efron', 
                  cluster = leaderID)
  
  results <- broom::tidy(model) %>%
    dplyr::select(estimate, robust.se)
  
  results 
}

cwe3_placebo_fe <- lapply(start_year, placebo_test_3_fe)

cwe3_placebo_fe_results <- bind_rows(cwe3_placebo_fe) %>%
  mutate(ll = estimate - (1.96*robust.se), 
         ul = estimate + (1.96*robust.se), 
         start_year: 1966:2006)


ggplot(cwe3_placebo_fe_results, aes(x = start_year, y = estimate, 
                                    ymin = ll, ymax = ul)) +
  geom_pointrange(size = .3) + 
  geom_hline(yintercept = 0, linetype = 'dashed') +
  scale_x_continuous(breaks = c(seq(1966, 2006, 1))) +
  scale_y_continuous(limits = c(-1.1, .85), 
                     breaks = c(seq(-1, .8, .2)), 
                     labels = c(seq(-1, .8, .2))) +
  labs(x = 'Window Start Year', y = 'Coefficient Estimate') +
  theme_pubclean() +
  theme(text = element_text(size = 12), 
        axis.text.x = element_text(angle = 75, vjust = .4),
        legend.position = 'bottom', 
        legend.key = element_blank())

ggsave('Fig5_CWE3_placebo_fe.pdf', 
       width = 9, height = 6)


# Supplemental Material Table A64 -----------------------------------------

stargazer(senior_cwe5_pcw, 
          se = list(senior_cwe5_pcw_se), 
          covariate.labels = c('CWE 5-year',
                               'Senior',
                               'Post-Cold War',
                               'ln(Oil/Gas Rents)',
                               'Personalism',
                               'Leader exit',
                               'Leader tenure',
                               'Leader tenure$^{2}$',
                               'Leader tenure$^{3}$',
                               'Election',
                               'Failed coup',
                               'ln(GDP per capita)',
                               'Cabinet size',
                               'Left truncated', 
                               'CWE 5-year X Senior', 
                               'Post-Cold War X Senior'),
          table.placement = '!hbt',
          label = 'tab:senior_cwe_pcw',
          dep.var.labels = 'Minister Exit',
          digits = 2,
          no.space = T,
          notes = 'Standard errors clustered by leader.', 
          star.cutoffs = c(0.05, 0.01, 0.001),
          title = 'Pooled Logit Model of Senior Minister Exit Across Cold War, Cold War End, and Post-Cold War from Figure 9',
          out = 'senior_cwe_pcw.tex')




